A mathematical model to serve as a clinical tool for assessing obstructive sleep apnea severity

Obstructive sleep apnea (OSA) is a sleep disorder caused by periodic airway obstructions and has been associated with numerous health consequences, which are thought to result from tissue hypoxia. However, challenges in the direct measurement of tissue-level oxygenation make it difficult to analyze the hypoxia exposure pattern in patients. Furthermore, current clinical practice relies on the apnea-hypopnea index (AHI) and pulse oximetry to assess OSA severity, both of which have limitations. To overcome this, we developed a clinically deployable mathematical model, which outputs tissue-level oxygenation. The model incorporates spatial pulmonary oxygen uptake, considers dissolved oxygen, and can use time-dependent patient inputs. It was applied to explore a series of breathing patterns that are clinically differentiated. Supporting previous studies, the result of this analysis indicated that the AHI is an unreliable indicator of hypoxia burden. As a proof of principle, polysomnography data from two patients was analyzed with this model. The model showed greater sensitivity to breathing in comparison with pulse oximetry and provided systemic venous oxygenation, which is absent from clinical measurements. In addition, the dissolved oxygen output was used to calculate hypoxia burden scores for each patient and compared to the clinical assessment, highlighting the importance of event length and cumulative impact of obstructions. Furthermore, an intra-patient statistical analysis was used to underscore the significance of closely occurring obstructive events and to highlight the utility of the model for quantitative data processing. Looking ahead, our model can be used with polysomnography data to predict hypoxic burden on the tissues and help guide patient treatment decisions.


DERIVATIONS AND FINITE DIFFERENCE FORMS OF EQUATIONS
Eqs. S12 (Inspiration) and S13 (Expiration) mmHg 104 mmHg (Guyton and Hall (2000)) P pa,O 2 Oxygen partial pressure in pulmonary arteries Eq. 3a mmHg Eq. 3a P pc,O 2 Oxygen partial pressure in pulmonary capillary compartment Eq. 3a mmHg Eq. 3a P pv,O 2 Oxygen partial pressure in pulmonary veins Eq. 3a mmHg 104 mmHg (Guyton and Hall (2000)) P sa,O 2 Oxygen partial pressure in systemic arteries Eq. 3a mmHg 104 mmHg (Guyton and Hall (2000) Guyton and Hall (2000)) cm 2 β p Oxygen solubility in alveolar-capillary membrane and blood plasma 1.4 × 10 −9 mol O 2 · mL −1 · mmHg −1 (Frank et al. (1997))       Figure S1. A flow diagram depicting the model physiology. The alveoli, pulmonary capillaries, and systemic capillaries were modeled as single compartments. Parameters related to the ambient air (external environment), anatomical dead space, heart, and body tissues were included in the model, but no mass balances were performed over these regions. All oxygen transfer in the systemic circulation was assumed to occur in the capillary region. Therefore, the pulmonary venous and arterial total oxygen concentrations were assumed to be equivalent to the systemic arterial and venous total oxygen concentrations, respectively.

Pulmonary Arterial Oxygenation
Equation 1, as restated below in its expanded form, was obtained by performing an unsteady-state mass balance over a moving control volume in the systemic capillary compartment, with an overall metabolic rate used for modeling oxygen transfer to the body tissues. To determine the mass transfer only occurring out of the control volume, the overall oxygen metabolic rate was multiplied by a ratio quantifying the fraction of the systemic capillary compartment occupied by the control volume: For the OSA patient cases, the systemic blood flow is a function of time as it varies with the input heart rate (Equation S1b ). In addition, since ∆t was made to vary proportionally with the heart rate, for data storage purposes (Appendix S5.1), ∆V remains constant. The finite difference form of Equation 1 is: Solving Equation S2 over time resulted in the following relation between the total oxygen concentration exiting and entering the systemic capillary compartment: Assuming that all oxygen transfer in the systemic circulation occurs in the capillary region, the total oxygen concentration in the pulmonary arteries and veins was related as: In addition, the total oxygen concentration in the pulmonary veins and that entering the systemic capillary compartment were related in time as: This assumes no oxygen leakage in the arteriole system. If there is possible arteriolar oxygen transfer to the tissues, the model will be adjusted accordingly. Due to their dependence on the systemic flow rate, t sys,cap , t cycle , and t sa were all functions of time for the patient cases. However, for the simulations, Q s , ∆t, t sys,cap , t cycle , and t sa were all constant due to an assumed constant heart rate (Table S3).  Figure S2. A schematic of the pulmonary capillary compartment. Deoxygenated blood enters from the pulmonary arteries and is oxygenated as a result of diffusion across the alveolar-capillary membrane. ∆z is the length of the differential control volume. Since the overall pulmonary flow passes through the compartment, an effective cross-sectional area is used (sum of all pulmonary capillaries), A eff , to determine the compartment length, L c . Oxygenated blood exits the compartment in the pulmonary veins.

Pulmonary Capillary Mass Transfer
Equation 4, as restated below, was obtained by performing an unsteady-state mass balance over a moving control volume in the pulmonary capillary compartment: S6b) For the OSA patient cases, the pulmonary blood flow is a function of time as it varies with the input heart rate (Equation S6b). Considering that the total pulmonary and systemic flow rates are equal to the cardiac output, ∆V for the pulmonary capillary compartment is constant due to the varying ∆t.
The finite difference form of Equation 4 is:  Figure S3. Effect of perfusion on model results. As the heart rate is increased, equilibration between the alveolar (dashed lines) and pulmonary capillary (solid lines) oxygen partial pressures occurs further along the capillary due to reduced residence time, indicating the limiting effect of perfusion in mass transfer.
To obtain an oxygen profile along the pulmonary capillaries, the starting boundary of the control volume was used to define its location. The pulmonary blood velocity was then used to track the movement of this volume along the pulmonary capillary compartment, introducing spatial dimensionality into the model and ensuring the significant influence of perfusion, which controls the equilibration between alveolar and end-capillary oxygen partial pressure (demonstrated with Figure S3).As a result, Equation S7 was transformed into: Due to its dependence on the pulmonary flow rate, v pc is also function of time for the patient cases. However, since ∆t varies proportionally with the heart rate, ∆z remains constant. Q p , ∆t, and v pc were kept constant for all the simulated cases due to an assumed constant heart rate (Table S3). Oxygen transfer into the control volume was modeled assuming a pseudo-steady state diffusion process across the alveolar-capillary membrane: In addition, the total oxygen concentration at the entrance of the pulmonary capillary compartment was set as pulmonary arterial value, while the total oxygen concentration in the pulmonary veins was set as the end-capillary value.

Alveolar Oxygenation
Converting Equation 5 to a finite difference form enabled efficient calculation of the alveolar oxygen partial pressure at any time point during inspiration. This equation was used under the condition that Converting Equation 6 to a finite difference form enabled efficient calculation of the alveolar oxygen partial pressure at any time point during expiration. This equation was used under the condition that ∆t was kept constant for all the simulations but varied for the OSA patient cases according to Equation S1c.

INITIAL CONDITION FOR PULMONARY CAPILLARY PROFILE
As an approximation, the initial oxygen profile (t=0) along the pulmonary capillary compartment was determined by solving a steady state mass balance over a fixed control volume.
The initial conditions for the total and dissolved oxygen concentrations in the pulmonary arteries were used for C pc,O 2 T (z = 0) and C pc,O 2 d (z = 0), respectively. The starting value for the total oxygen in the pulmonary arteries was calculated using Equation S4. As an additional form of validation, the predicted average oxygen partial pressure in the pulmonary arteries (P pa,O 2 = 41 mmHg) approximates the expected physiological value of 40 mmHg (Guyton and Hall (2000)) for a normal subject at rest.

ASSUMPTION OF PSEUDO-STEADY STATE DIFFUSION ACROSS ALVEOLAR-CAPILLARY MEMBRANE
The molar flux of oxygen across the alveolar-capillary membrane can be represented as (Popel (1989)): was used as an approximation of the effective diffusion coefficient, which was applied to estimate the diffusion time across the alveolar-capillary membrane: The process time (time for boundary changes) was taken to be on the order of the time for one breath. With a normal breathing rate of 12 breaths/min, this would be 5 s, while, with a hyperventilatory rate of 24 breaths/min, the process time would be 2.5 s. For both of these cases, the process time was concluded to be sufficiently larger than the diffusion time, allowing for the pseudo-steady state simplification.  Figure S4. Patient 1 raw and normalized nasal pressure signals over isolated normal breathing region.

ESTIMATION OF OSA PATIENT LUNG VOLUME USING RECORDED NASAL PRESSURE
For each patient, the recorded pressure signal was observed to have a non-zero mean, which resulted in unphysiological values for the lung volume after conversion. Therefore, to allow for a more realistic representation of the volume, each pressure signal was normalized to a mean of zero by subtracting a moving average taken over 70 seconds for Patient 1 and 50 seconds for Patient 2. The time for each moving average was chosen such that it was on the order of a few breaths, while ensuring that the normalized signal maintained the integrity of the raw pressure signal (Figures S4 and S5).
The following equations were used to simulate normal breathing patterns and corresponding flow rates for the OSA patients, with the tidal volume and breathing rate for each determined using the patient  Figure S6. Actual and simulated signals for Patient 1. A) Normalized patient nasal pressure signal during isolated normal breathing sequence. B) Simulated normal lung volume and flow rate using breathing rate from isolated patient signal (b r = 20 breaths/min) and tidal volume from patient height (V T = 0.51 L).
height-based ideal body weight and normalized pressure signal, respectively: and  Figure S7. Actual and simulated signals for Patient 2. A) Normalized patient nasal pressure signal during isolated normal breathing sequence. B) Simulated normal lung volume and flow rate using breathing rate from isolated patient signal (b r = 12 breaths/min) and tidal volume from patient height (V T = 0.56 L).
A positive nasal pressure indicated inspiration, while a negative nasal pressure indicated expiration. Therefore, for each patient, the fitting parameter for conversion to nasal flow was determined by setting the average maximum nasal pressure (after normalization) during the isolated normal breathing sequence ( Figures S6A and S7A) to the maximum simulated inspiratory flow rate ( Figures S6B and S7B): Each normalized nasal pressure signal was converted to nasal flow as: The patient alveolar volume, approximated as the lung volume, was determined by integration of the obtained flow signal:  Figure S8. Pressure swings during a segment with obstructive events. A) Raw pressure signal during polysomnography study. B) Simulated lung volume after converting normalized pressure using fitting parameter. The sections outlined in red show a sequence with multiple obstructive events and drastic pressure swings, resulting in a low lung volume after conversion. The grayed portions indicate obstructive respiratory events, as labeled in (A).
While analyzing Patient 2, a region with multiple obstructive events and drastic pressure swings was observed ( Figure S8A, outlined in red). Although the majority of the volume signal after conversion is reasonable, this section is an example where relying solely on nasal pressure results in values of the simulated lung volume that are much lower than the observed residual volume in OSA patients (Abdeyrim et al. (2015)), with the lowest point dropping close to 1.0 L ( Figure S8B). Using RIP chest and abdomen signals, which are available within most clinical polysomnographies, in conjunction with the nasal pressure will likely provide a more accurate representation of patient lung volume in such cases.

Pulmonary Arterial Oxygenation
Output for next time step Figure S9. Model framework. Breathing and heart rate data (patient cases), physiological and solving parameters, and initial conditions were input to the model. Each equation was run for multiple time iterations to obtain the outputs.
All MATLAB codes used to solve the model equations can be found at GitHub: https://github.com/Cardiovascular-Modeling-Laboratory/OSAModel.git.

Equation S4
was implemented into the MATLAB code as: Here, the for loop index, i, was used to represent a single time iteration. For the patient cases, dt was made to vary proportionally with the heart rate to maintain the number of iterations defining the movement of the control volume through the entire systemic circulation, artery/arteriole system, and systemic capillaries. This was done to avoid any gaps in the output data arrays for C pa,O 2 T , C sa,O 2 T , and C sv,O 2 T , in the case of a decrease in heart rate between subsequent time iterations, or overwriting of previous data, in the case of an increase in heart rate. In Equation S22a, dt(1) and t cycle (1) represent the initial conditions for dt and t cycle , respectively. The floor t cycle (i) dt(i) term in Equation S22 was used to determine the number of iterations defining the movement of the control volume through systemic circulation. Due to the variation of dt for the patient cases, this value remains constant. Similarly, Equations S5 and S3 were implemented as: The floor t sa (i) dt(i) and floor t sys,cap (i) dt(i) terms were used to determine the number of iterations defining the movement of the control volume through systemic artery/arteriole system and systemic capillaries, respectively. These values also remain constant due to the varying dt in the patient cases. For all the simulations, Q s , dt, t sys,cap , t cycle , and t sa were implemented as constants in the code.

Pulmonary Capillary Mass Transfer
Equation S8 was implemented into the MATLAB code using the following setup: The result of C pc,O 2 T was a matrix, where each column represented the spatial oxygen profile at a specific time point. In addition, the total oxygen concentration in the pulmonary veins was determined as: For all the simulations, Q p was implemented as a constant in the code.

Alveolar Mass Transfer
Equations S12 and S13 for inspiration and expiration, respectively, were implemented into the MATLAB code as: and These equations were coded into a separate function. For all the simulations, dt was implemented as a constant in the code, but, for the OSA patient cases, it was varied according to Equation S22a.

Time Step for Simulations and Patient Cases
The initial time step (dt(1)) for the simulations and patient cases was chosen based on the following conditions: ensuring i) model stability over the period of study, ii) model convergence, and iii) maximal computational efficiency. Convergence was tested using the normal breathing simulation. The model was run at values of the time step lower than the selected value of 0.01 s (Table S3), and it was ensured that the solutions were reasonably close for each one. dt(1)>0.01 s was not chosen as the solution began to lose stability and diverge. Although the model can be run at dt(1)<0.01 s, computation time is significantly longer. Therefore, 0.01 s was chosen to maintain a balance between all three conditions. For the patient cases, since dt(1) varies with the heart rate (Equation S22a), the initial time steps of 0.008 s for Patient 1 and 0.007 s for Patient 2 were chosen such that the maximum value would be 0.01 s:

Estimation of OSA Patient Lung Volume Using Recorded Nasal Pressure
For Patient 1, the raw nasal pressure signal was normalized as: For Patient 2, the raw nasal pressure signal was normalized as: dt m was the time step for measurement. Eq. S21 was implemented into the MATLAB code for each case to determine the time-dependent lung volume using cumulative integration: Start i and End i represented the starting and ending time elements for each breathing pattern, respectively, and were chosen to avoid the time period when CPAP was administered to the patients during the sleep studies as well as segments of time with any signal loss in recorded data due to patient movement. In addition, Start i was chosen around a point where the nasal pressure was zero, before inspiration, which followed an ideal normal breath with similar inspiratory and expiratory flow rates. The resulting vector of V A,OSA for each patient was imported into the main code to serve as a column for a reference lookup table.

Lookup Tables
Interpolation using reference lookup tables was employed to determine the values of alveolar volume, heart rate (for patient data), and dissolved oxygen concentration at each time step. For the OSA simulations, variations of Equation 7 were used to construct breathing patterns. The time and alveolar volume data for each pattern was converted into a .mat file and imported into the main code. Heart rate was kept constant for the simulated cases. A similar conversion process was repeated for the patient data, which also included the variable heart rate input in the .mat file. The data extracted from the .mat files served as reference lookup tables in the main code. Interpolation for alveolar volume and heart rate was carried out using: and HR(i) = interp1(t data , HR data , t(i)) (S36) t data (time) represented the input column of the lookup table, while V data (alveolar volume) and HR data (heart rate) were the output columns.
Due to the non-linearity of Equation 3 and to minimize computational effort, a lookup table was constructed to determine the dissolved oxygen concentration at each spatial position. Equation 3 was used to calculate the hemoglobin oxygen saturation for a range of dissolved oxygen concentrations. To ensure accuracy, the dissolved oxygen column in the table ranged from 0 to 2×10 −6 mol/mL with increments of 2.5×10 −10 mol/mL. These values were then used with Equation 2 to determine the corresponding total oxygen concentrations. This allowed the use of interpolation: CT (total oxygen concentration) represented the input column of the lookup table, while Cd (dissolved oxygen concentration) was the output column. j in Equation S37 defines the spatial location. j = pa, pc, pv, sa, and sv for the pulmonary arteries, pulmonary capillary compartment, pulmonary veins, systemic arteries, and systemic veins, respectively.

Comparison of Normal Breathing Results with Expected Physiological Values
Expected physiological values for the dissolved oxygen concentration in the blood vessels are not commonly found in literature. Therefore, after converting the wide range of blood oxygen partial pressures in the systemic vessels (Guyton and Hall (2000); Ortiz-Prado et al. (2019);van Faassen et al. (2009)) using β p (Table S2), values of 119-140 µM and 38-63 µM were obtained for the systemic arteries and veins, respectively. The average values predicted by the model (Figure 1) for the systemic arteries and veins fall within these ranges. It is, however, important to note that the dissolved oxygen is sensitive to the choice of solubility factor. Percent decreases in saturation and dissolved oxygen concentration for each simulation were calculated by comparing the average normal values to the absolute minimum points over the time of study:

Calculation of Percent Decreases for OSA Simulations
t OSA is the duration of the simulated OSA breathing pattern, while t total is the total time for breathing data in each case, both in seconds.

Calculation of Proposed Burden Scores
In Equation S41, n was defined as the total number of wake sequences considered (n = 3 for the current analysis). The subscripts s and e in Equation S41a represent the start and end times for each wake sequence, respectively. The proposed hypoxia burden score was determined as: In Equation S42, t 1 represents the start time for the analyzed sequence (t 1 = 360 s for the total sequence), and t 2 represents the end time for the analyzed sequence. For the proposed burden score of the hypoxic period over the total patient sequence, the upper threshold was calculated as: C sa,O 2 d{std,wake} represents the standard deviation of the arterial dissolved oxygen during wakefulness.  Figure S10. Schematic of area used to calculate proposed burden score for the hypoxic periods over the total sequences.
The proposed burden score for hypoxic period during the total sequence was calculated using the area deviation of the dissolved arterial oxygen from C sa,O 2 d{hypoxia} ( Figure S10). To determine the score, this area value was normalized by the total time for hypoxia. To evaluate the clinical applicability of the model over longer time periods, data from 3-hour and 2-hour portions of baseline polysomnographies, without any CPAP administration, was fed to the model for Patient 1 and 2, respectively. The heart rate and lung volume were used as direct inputs (Figures S11 and S13), while hemoglobin oxygen saturation and dissolved oxygen concentration in the systemic arteries and veins were output from the model (Figures S12 and S14).  The differences for the one-sample t-tests were calculated by subtracting the overall average of each variable from averages of the same variable over identified event series intervals: